set.seed(1)
#if (!requireNamespace("BiocManager", quietly = TRUE))   
#  install.packages("BiocManager",repos = "http://cran.us.r-project.org")  
#BiocManager::install("fgsea")
#BiocManager::install("limma")
#BiocManager::install("Biobase")
#install.packages("ggplot2")
#install.packages("igraph")
library(limma)
library(Biobase)
library(ggplot2)
library(igraph)
library(data.table)
library(fgsea)
library(readr)
library(biomaRt)
library(edgeR)

##############################
# Number of mRNA
n_mRNA = 19504

bonobo = data.frame(fread("/home/esaha/BONOBO/network/miRNA_mRNA_breastCancer/indegree_mRNA.txt"))
genes = bonobo$V1

indegree = data.frame(bonobo[-(1:n_mRNA),-1])
rownames(indegree) = genes[-(1:n_mRNA)]

# Get phenotypic data
phenotype = read.csv("~/BONOBO/data/miRNA_mRNA_breastCancer/GSE19783_phenotype.txt", sep="")

subtypes = phenotype$breast.cancer.subtype.ch1[match(colnames(indegree), rownames(phenotype))]
subtypes[is.na(subtypes)] = "not classified"
#subtypes[subtypes == "not classified"] = "unknown"

# limma model to identify difference between each pair of subtypes
subtypes = factor(subtypes)
# Table for how many significant genes are differentially coexpressed between pair of subtypes
table_sigGene = matrix(0, length(levels(subtypes)), length(levels(subtypes)))
rownames(table_sigGene) = gsub("subtypes", "", levels(subtypes))
colnames(table_sigGene) = gsub("subtypes", "", levels(subtypes))
for (i in 1:length(levels(subtypes))){
  subtypes0 = relevel(subtypes, ref = levels(subtypes)[i])
  design = model.matrix(~ subtypes0)
  fit <- lmFit(indegree, design)
  fit <- eBayes(fit)
  for (j in 1:ncol(fit)){
    if (j != i){
      coef_name = paste("subtypes0", levels(subtypes)[j], sep = "")
      tb = topTable(fit,coef=coef_name,number=Inf)
      table_sigGene[i,j] = length(tb$P.Value[which(tb$P.Value < 0.05)])
    }
  }
}
rownames(table_sigGene) = levels(subtypes)
colnames(table_sigGene) = levels(subtypes)

# Make heatmap comparing all groups

library(ggplot2)
library(gplots)
library(RColorBrewer)

mycol2 <- colorRampPalette(rev(brewer.pal(11,"RdBu")))(50)
pdf("/home/esaha/BONOBO/plots/mRNA_miRNA_breastCancer/GSE19783_subtype_difference_miRNAmRNA_correlation.pdf",h=25,w=16)
heatmap.2(as.matrix(-table_sigGene),density.info="none",trace="none",col=mycol2,symbreaks=T,symkey=T, key = F,
          cexRow=2.5, cexCol=2.5, srtCol = 45, keysize=0.5, margins = c(20, 15), key.title=NULL, key.xlab="NES", cellnote=table_sigGene, notecol="black", notecex=2)
dev.off()
